# File: MV.R
# R source code for Markowitz Protfolio optimization.
# Author: Congxing Cai (congxing@stanford.edu)
optimal.weights <- function(mu, sigma, aversion) {
	# Compute the optimal allocation weights for the next period.
	# Args:
	#   mu: N*1 vector, of the expected returns of the portfolios
	#   sigma: N*N matrix, of the covariance of the portfolios
	# if (mu.expected == 0) {
		# w <- weights.mvp(mu, sigma)
	# } else {
		# w <- optimal.plugin(mu, sigma, as.numeric(mu.expected))
	# }
	w = (1/aversion) * solve(sigma) %*% mu
	return(w)  
}

minvar.mvp <- function(mu, sigma) {
	# Compute the minimum-variance point using plugin matrix: equation (3.13).
	# Args:
	#   mu: N*1 vector, of the expected returns of the portfolios
	#   sigma: N*N matrix, of the covariance of the portfolios
	iSigma <- solve(sigma)
	ind <- as.matrix(rep(1, length(mu)))
	A = as.numeric(t(ind) %*% iSigma %*% mu)
	C = as.numeric(t(ind) %*% iSigma %*% ind)
	return(A/C)
}

weights.mvp <- function(mu, sigma) {
	# Compute the minimum-variance weight vector using plugin matrix: equation (3.14)
	# Args:
	#   mu: N*1 vector, of the expected returns of the portfolios
	#   sigma: N*N matrix, of the covariance of the portfolios
	iSigma <- solve(sigma)
	ind <- as.matrix(rep(1, length(mu)))
	C = as.numeric(t(ind) %*% iSigma %*% ind)
	return(iSigma %*% ind / C)
}


optimal.plugin <- function(mu, sigma, mu.expected) {
	# Compute Markowitz optimal weights, given expected return, using plugin matrix: equation (3.11) - allowing short selling
	# Args:
	#   mu: N*1 vector, of the expected returns of the portfolios
	#   sigma: N*N matrix, of the covariance of the portfolios
	#   mu.expected: the expected return
	iSigma <- solve(sigma)
	ind <- as.matrix(rep(1, length(mu)))

	# Lagrange multipliers under short selling
	B = as.numeric(t(mu) %*% iSigma %*% mu)
	A = as.numeric(t(ind) %*% iSigma %*% mu)
	C = as.numeric(t(ind) %*% iSigma %*% ind)
	D = as.numeric(B %*% C - A %*% A) 

	w = (B * iSigma %*% ind - A * iSigma %*% mu + mu.expected * (C * iSigma %*% mu - A * iSigma %*% ind)) / D
	return(w)
}


get.efficient.frontier <- function(mu, sigma, mu.targets) {
	# Compute the (sigma, mu) pairs for the mu.targets to construct the efficient frontier
	# Args:
	#  	mu: N*1 vector, of the expected returns of the portfolios
	#   sigma: N*N matrix, of the covariance of the portfolios
	#   mu.targets: the samples of targeted returns.   
	x <- numeric(length(mu.targets))
	y <- numeric(length(mu.targets))
	for (i in 1:length(mu.targets)) {
		mu.target <- mu.targets[i]
		w <- optimal.plugin(mu, sigma, mu.target)
		y[i] <-	t(w) %*% mu
		x[i] <- sqrt(t(w) %*% sigma %*% w) 	
	}	
	return(data.frame(x, y))
}

#### Experiments ####

# For section 1: Introduction
mv.plugin.intro <- function(){	
	# Markowitz M-V optimization using plugin estimates of mu and sigma.
	ret <- data.frame(FF6$smlo_vwret, FF6$smme_vwret, FF6$smhi_vwret, FF6$bilo_vwret, FF6$bime_vwret, FF6$bihi_vwret)
	names(ret) <- c("smlo", "smme", "smhi", "bilo", "bime", "bihi")
	mu <- as.matrix(mean(ret))
	sigma <- cov(ret)
	w.mvp <- weights.mvp(mu, sigma)
	png("../results/markowitz_enigma.png")
	barplot(as.vector(w.mvp), main="Portfolio Weights using Plug-in Estimates", xlab="portfolio", ylab="weight", names.arg=c("smlo", "smme", "smhi", "bilo", "bime", "bihi"), col=c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00"), legend = c("smlo", "smme", "smhi", "bilo", "bime", "bihi"))
	dev.off()
	
	cap.wt <- get.market.capital.weight(FF6)
	names(cap.wt) <- c("smlo", "smme", "smhi", "bilo", "bime", "bihi")
	cap.wt <- as.matrix(cap.wt[1,])
	png("../results/cap_weights.png")
	barplot(as.vector(cap.wt), main="Market Capital Weights of the Portfolios", xlab="portfolio", ylab="weight", names.arg=c("smlo", "smme", "smhi", "bilo", "bime", "bihi"), col=c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00"), legend = c("smlo", "smme", "smhi", "bilo", "bime", "bihi"))  
	dev.off()
}
